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Abstract 

We review a family of local algorithms that permit the simulation of charged particles with purely local dynamics. 
Molecular dynamics formulations lead to discretizations similar to those of "particle in cell" methods in plasma 
physics. We show how to formulate a local Monte-Carlo algorithm in the presence of the long ranged Coulomb 
interaction. We compare the efficiencies of our molecular dynamics and Monte-Carlo codes. 
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1. Introduction 

Computer modeling of charged systems is de- 
manding due to the range of the Coulomb inter- 
action [1]. The direct evaluation of the Coulomb 
sum for N particles, U c = J2i<j e i e j I '^ e o r ij ) re- 
quires computation of the separations ry between 
all pairs of particles, which implies 0(N 2 ) opera- 
tions are needed per sweep or time step. Most large 
scale codes in use at the moment find the elec- 
trostatic energy and forces using the fast Fourier 
transform [2] after interpolating charges to a mesh. 
For an ensemble of charges interpolated to O(M) 
nodes of a lattice this gives the electrostatic en- 
ergy in a time which scales as 0(N + M In M) on 
a mono-processor machine. 

It is rather surprising that methods that are 
based on solution of the Poisson equation are so 
dominant in mesh based molecular dynamics. In 
plasma physics, one may use a global solver or al- 



ternatively one integrates the full local Maxwell 
equations for the electromagnetic field. We shall 
show in this article that a similar approach based 
on Maxwell's equations also gives good results in 
molecular dynamics; Coulomb's law is implied by 
Maxwell's equations in the quasi-static limit. The 
advantage of this approach is that parallel imple- 
mentations which exhibit 0(N + M) scaling are 
rather easy to write. 

After abstracting the part of Maxwell's equa- 
tions which is essential for Coulomb's law we intro- 
duce a local Monte-Carlo algorithm based on an 
auxiliary electric field E. This method also requires 
an effort of 0(N + M) per sweep and leads to sub- 
stantial speedups compared to conventional molec- 
ular dynamics methods. Again the algorithm is eas- 
ily parallelizable. Introduction of cluster moves for 
the field E give further increases in efficiency in the 
O(M) mesh work; this is most useful when N/M 
is small. 
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2. Particle in Cell codes in Plasma Physics 

Particle in cell methods [3,4] are widely used in 
the simulation of plasmas coupled to Maxwell's 
equations. Applications are widespread and in- 
clude astrophysics, laser physics and fusion re- 
search. In these algorithms the continuum mag- 
netic and electric fields are discretized on mesh. 
Particles, however, still move in the continuum:- 
just like mesh based molecular dynamics algo- 
rithms. The principle technical difficulty lies in 
correctly coupling the field and particle degrees 
of freedom. The coupling is usually performed by 
interpolation of both the electric charges and the 
electric currents to the mesh and then by a back 
interpolating the electrodynamic forces from the 
mesh onto the particles. 

Typically the charge of the particles, d is inter- 
polated to the nodes while the current, Ji = is 
assigned to the links. This double interpolation is 
difficult; conservations laws may not be preserved. 
The importance of exact implementation of conser- 
vation laws is now widely recognized from the field 
of geometric integration [5]. Many charge/current 
interpolation schemes do not satisfy the continuity 
equation 

^+divJ = (I) 

We shall show below that the breakdown of con- 
tinuity leads to errors in the electrostatic interac- 
tions between charges. 

Until now the principle solutions to this prob- 
lem of charge continuity were either the use of a 
low order (and thus noisy) interpolation scheme for 
which continuity is exact [3] or the use of higher 
order schemes with frequent global corrections in 
which the longitudinal field components are regu- 
larly thrown out and then reinitialized from solu- 
tion of the Poisson equation. A recent and inter- 
esting hybrid approach considers modifications to 
Maxwell's equations [6,7] which stabilize their so- 
lutions even in the presence of inconsistencies in 
the interpolation scheme. 

The use of a lattice to interpolate the field de- 
grees of freedom leads to two other sources of error 
in the energy of a system of charged particles in the 
quasi-static limit. Firstly the pair interaction be- 



tween two particles is given by a lattice Green func- 
tion, rather than the exact interaction in 1/r. In 
the simplest discretizations this leads to an inter- 
action potential of the form 1/r + 0(a 2 /r 3 ) where 
a is the lattice spacing. A second and more seri- 
ous artefact is the variation of the self-energy of 
the particles as a function of their position with 
respect to the lattice. This leads to a periodic 1- 
body potential for a particle i with an amplitude 
which varies as ef/eoa [8]: As the lattice spacing 
decreases this energy variation diverges so the dis- 
cretization becomes worse. In the plasma literature 
this 1-body energy is often ignored, one simply uses 
a very coarse mesh with large a so that the kinetic 
energy is large compared with this energy scale; in 
a thermalized plasma the mesh size must be larger 
than the Bjerrum length, If, — e 2 / 'A-KeksT . 

In mesh based molecular dynamics codes these 
two artefacts in the energy are also present. Usually 
the errors are reduced by convolving the charge of 
each particle over a large number of lattice points 
[2,9,10], or equivalently by multiplying the struc- 
ture factor in Fourier space by an appropriate shap- 
ing factor. We need to implement similar meth- 
ods if we wish to use particle in cell techniques for 
molecular dynamics. 

2.1. Implementation for Molecular Dynamics 

Recently two groups have implemented a molec- 
ular dynamics algorithm based on particle in cell 
discretizations. As noted above conventional for- 
mulations arc not quite good enough for direct ap- 
plication in molecular dynamics. This new work 
has 

- Shown how Maxwell's equations give Coulomb's 
law 

- Led to new charge conserving integrators 

- Controled lattice artefacts in the energy 

In our own work [11] we used interpolation al- 
gorithms based on polynomial splines in order to 
'tune' the error in the 1-body and pair potentials 
in a manner very similar to that used in conven- 
tional molecular dynamics codes based on Poisson 
solvers [9]. To do this we generalized the charge 
conserving algorithm of [3], valid for linear inter- 
polation to general B-spline interpolation. The re- 
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suiting algorithm [11] conserves charge to machine 
accuracy for finite time steps and is time reversible. 
In order to simplify the implementation we chose 
to drop the Lorentz force from the particle equa- 
tions of motion which leads to non-Hamilton (and 
thus non symplectic) dynamics; however the dy- 
namics do still conserve the total volume in the ap- 
propriately defined configuration space. This con- 
servation law implies the existence of a Gibbs dis- 
tribution for the combined field-particle system. 
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Fig. 1. Cpu time for the field integration (□), field-particle 
couplings (A) and total time (()) in a system of size L = 30a 
as a function of N on an single AMD Athlon CPU. Also 
shown arc results for the same system treated with the 
P3M method using a charge interpolation of the same order 
(filled symbols), ■ measures the total time in Fourier based 
work. 

Our molecular dynamics code has been run in 
two distinct modes. Firstly a mode in which the 
transverse components of the electric field and 
the particles are coupled to a thermostat. For this 
mode we are able to demonstrate, analytically, 
convergence to a Boltzmann- Gibbs distribution 
of particles interacting as if with a instantaneous, 
non-retarded Coulombic interaction. In the second 
mode we continue to thermostat the particles but 
anneal the transverse components of the electric 
field to zero temperature. This mode is very simi- 
lar to Car-Parrinello [12] simulations in quantum 
molecular dynamics. Here it is the ratio of the 
typical particle velocity, v to the speed of light, c 
which is tuned as the optimization parameter. We 
expect (but here have no proof) that the correct 
Gibbs distribution is generated when v/c is suffi- 
ciently small. We again verified numerically that 
the correct effective interaction is found. A com- 



parison of the speed of our code compared with a 
conventional Fourier based Poisson solver [13] is 
given in Fig. (1). At very low densities our code has 
a slight advantage, at higher densities the conven- 
tional algorithm behaves better. The advantages 
in each regime however are rather modest. 

Pasichnyk and Duenweg [14] have taken a rather 
different and more direct route to eliminating the 
errors in the potentials. They make the remark 
that for low order interpolation schemes one can 
rather easily calculate the exact numerical value of 
the 1-body potential. They then perform a direct, 
numerical subtraction of this energy in their code. 
They have performed a detailed comparison of this 
simple and elegant solution with simulations per- 
formed with conventional Fourier methods. They 
find that their code is competitive in both speed 
and accuracy with conventional methods. The au- 
thors work again in the limit v/c small but with 
no thermostating of the electric field. In contrast 
with our own work they have a full parallel imple- 
mentation using MPI. 



3. Origin of Coulomb's law 

In electromagnetism Coulomb's law comes 
[15] from a local expression for the energy: U = 
J £oE 2 /2 d 3 r, and the imposition of Gauss' law 
div E — p/e = 0. This motivates the use of the 
following partition function for the electric field 



Z(p) = J DV 11 5 (div E - p/eo) e 



U/k B T 



(2) 



We change integration variables to e = E + V</> 
with V 2 </> = —p/eo an d find 

Z = jDeS(drve)e-^I^- V ^ d3r . (3) 

The cross term in the exponential is zero by inte- 
gration by parts: / V</> • e d 3 r = 0, so that 

z = e -^j m f J DeS{dive)e -^fe 2 * 3 r (4) 

We conclude that Z = Zcoui x const, with ZqouI 
the partition function of particles interacting with 
the conventional Coulomb expression: Relative 
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statistical weights are unchanged if we allow the 
transverse field to vary freely, rather than quench- 
ing it to zero. To sample the partition function 
Eq. (2) we have to find solutions to Gauss' law. 
Integrating Maxwell's equations is one method of 
continuously and locally updating the system in 
such a way that Gauss' law is always satisfied. 

We now understand the importance of charge 
conservation in the generation of Coulomb's law: 
Consider the Maxwell equation 

e — = -J + curlH (5) 

and take its divergence, div (e E + J) =0. If we 
impose continuity Eq. (1) then we conclude that 

^(eodivE-p) = (6) 

and Gauss' law is valid. In the absence of continu- 
ity Gauss' law is violated and from the above argu- 
ment Coulomb's law is not generated as the effec- 
tive thermodynamic interaction between charges. 



4. Monte-Carlo sampling 

As well as performing molecular dynamic sam- 
pling of Eq. (2) we have performed Monte-Carlo 
sampling of the constrained partition function for 
the case a simple lattice gas [16] and also [8] for off- 
lattice models. These codes work with a mixture of 
two different Monte-Carlo moves. The first is the 
motion of a particle together with a slaved update 
of the electric field in the vicinity of the particle. 
The second is a grouped update of the electric field 
on the four links defining a plaquette of the mesh. 
The two moves are randomly mixed so that the 
probability of plaquette updates is p. These two 
moves have been related to the two terms on the 
right hand side of Eq. (5) [17]. 

We now compare molecular dynamics and 
Monte-Carlo codes for efficiency: a priori Monte- 
Carlo should be faster. 

- One does not restrict step sizes in order to ensure 
accuracy and stability of the dynamics. 

- Energy is easier to calculate than force, partic- 
ularly when interpolating to and from a mesh. 
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Fig. 2. Ratio of single particle diffusion coefficients in a sim- 
ple electrolyte obtained from molecular dynamics Dmd, 
and Monte Carlo -Dmc; time is measured in sweeps for 
each code. The particles also interact with a truncated, re- 
pulsive Lennard Jones potential of range r c = 2 1 / 6 <r, where 
<7 = 2a. Measurements as a function of relative interac- 
tion strength for several volume fractions: <3? = 0.15 (I), 
<S> = 0.3 (♦), * = 0.45 (A), * = 0.6 (★). 

We explore the first point in Fig. (2) by comparing 
the single particle diffusion coefficients in a sim- 
ple electrolyte as a function of interaction strength 
(Bjerrum length) using molecular dynamics pa- 
rameters from [10]. We also added a truncated, re- 
pulsive Lennard Jones potential to prevent over- 
lap of the particles. Simulation time is measured 
in sweeps, and the interaction strength is varied 
by changing the Bjerrum length lb from a/2 (weak 
Coulomb interaction) to 5tr (strong interaction). 
We find that in both codes the diffusivity decreases 
with increasing but the ratio remains roughly 
constant. At low volume fractions the Monte Carlo 
code is up to 80 times faster than the molecular 
dynamics code. This advantage begins to weaken 
at higher volume fractions, when collisions (due to 
the LJ potential) between particles slow down the 
diffusion apparently more strongly in MC. Never- 
theless Monte Carlo always retains an advantage. 
We have not evaluated the second factor in the rel- 
ative efficiencies since our two codes are structured 
very differently. 

Until now we have only considered media with 
uniform dielectric properties however Maxwell's 
equations are valid for arbitrary spatial variation 
of e(r). They remain local and thus give an O(N) 
algorithm for treating implicit solvents in a situa- 
tion where conventional Fourier codes break down. 
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The Monte-Carlo algorithm can also be applied to 
such media [18] . To do this one works with the elec- 
tric displacement D rather than the electric field. 
The Gauss law constraint now reads div D = p, 
the local expression for the energy is given by U = 
/D 2 /2e(r)d 3 r. 

4.1. Effective Monte-Carlo Field Dynamics 

Our Monte-Carlo algorithm is based on local up- 
dates to both the particle positions and the elec- 
tric field. Both the particles and the electric field 
exhibit diffusive dynamics, in contrast to the prop- 
agative behaviour of free particles and waves in the 
Maxwell equations. We have shown that the field 
dynamics are given by the Langevin equation [17]: 

-q£=D (V 2 E - grad p/e ) - 3/e + £(t) (7) 

which replace Maxwell's equations in our Monte- 
Carlo scheme. Eq. (7) is particularly interesting 
when there are free, mobile charges so that the 
electric current is linked to the electric field via the 
equation J = crE. In this case the dispersion law of 
the transverse electric field develops a gap so that 
iuj = cr/eo + Dq 2 

The relative diffusion rates of the electric and 
charge degrees of freedom can be adjusted by mod- 
ifying p and thus D. We initially believed that the 
algorithm should be run in a regime in which the 
effective diffusion coefficient of the electric field, D 
is comparable to or somewhat larger to that of the 
particles: The particles are then always interact- 
ing via a field that is close to equilibrium and are 
not slowed down by the "drag" from the electric 
degrees of freedom. 

Numerical experimentation has shown that one 
should use a much lower rate of update. There are 
several ways of understanding this surprising effi- 
ciency of the code 

- The gap in the electric dispersion law gives a 
fast relaxation of the electric degrees of freedom 
even at very long wavelengths. 

- In certain limits Coulomb interactions are gen- 
erated even if the plaquette updates are sup- 
pressed entirely (p = 0): a most surprising and 
deep result linked with the statistical mechanics 
of discrete plaquette models [19]. 



4.2. Cluster Monte- Carlo for E 

The efficiency of the Monte-Carlo algorithm is 
in part due to the gap which occurs in the disper- 
sion law for the transverse electric field. This gap 
has as its origin the finite conductivity of a system 
with a finite density of free charges. What happens 
when this density is small or zero? Are we obliged 
to spend much more of our computer budget up- 
dating the plaquettes in order to increase D and 
thus sample the O(M) electric degrees of freedom? 
For N/M small this would reduce the effective effi- 
ciency of the algorithm. Recently we have found a 
way [20] of reintroducing a gap in the electric field 
spectrum even when a = based on ideas used to 
accelerate the simulation of quantum spin models 
[19]. This quantum algorithm is in turn inspired by 
the Swendsen-Wang algorithm for the Ising model. 
In the language of electrostatics this new cluster 
algorithm works by nucleating a pair of virtual par- 
ticles of opposite sign. These particles then evolve 
by kinetic Monte-Carlo until they annihilate tak- 
ing the system back to a state of no free parti- 
cles. The virtual particles re-establish the gap in 
the electric spectrum and lead to resampling of the 
electric field in O(l) sweeps. The algorithm gives a 
dynamic exponent for the field z — 0, rather than 
z = 2 which is implied by Eq. (7) when a = 0. 



5. Perspectives 

We have introduced a series of algorithms for 
the local simulation of particles interacting via 
electrostatic interactions. Molecular dynamics im- 
plementations are broadly similar in performance 
to conventional Fourier based solvers. They can, 
however, be trivially parallelized in an environ- 
ment with limited intcrprocessor bandwidth; this 
is much harder to do with Fourier based solvers. 

We have also formulated an off-lattice Monte- 
Carlo algorithm that is faster than the molecular 
dynamics equivalent while keeping the advantage 
of locality. Finally analogies with "worm" or quan- 
tum cluster algorithms allow one to integrate over 
the field degrees of freedom with remarkable effi- 
ciency. 
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We note that the above discretizations of elec- 
trostatics are very similar to those introduced in 
the 19'th century. At this time there was a very 
poor intuition as to the physical content and im- 
plications of Maxwell's equations. FitzGerald in- 
troduced a series of mechanical models of the ether 
which are essentially conserving spatial discretiza- 
tions of the continuum Maxwell equations. In these 
models a series of wheels are connected by elastic 
bands [21]. The rotational motion of the wheels is 
equivalent to the magnetic field while the exten- 
sion, E of the elastic bands is equivalent to the elec- 
tric field. If we perform Monte-Carlo on FitzGer- 
ald's wheels we again recover our local algorithm 
for Monte-Carlo simulation. 
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